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Using a microfluidics device filled with a colloidal suspension of microspheres, we test the laws of 
diffusion in the limit of small particle numbers. Our focus is not just on average properties such as 
the mean flux, but rather on the features of the entire distribution of allowed microscopic trajectories 
that are possible during diffusive dynamics. The experiments show that: (1) the flux distribution is 
Gaussian; (2) Fick's Law — that the average flux is proportional to the particle gradient — holds 
even for particle gradients down to one or zero particles; (3) the variance in the flux is proportional 
to the sum of the particle numbers; and (4) there are backwards flows, where particles flow up a 
concentration gradient, rather than down it, and their numbers are well-predicted by theory and 
consistent with a new Flux Fluctuation Theorem. 



I. INTRODUCTION 

Fick's Law, which describes the diffusion of atoms, molecules, and particles, is important in many areas of science, 
and is the basis for engineering models of material transport. One statement of Fick's First Law is that the average 
particle flux is proportional to the average concentration gradient , 

where (J) is the observed macroscopic flux and (c) is the concentration of particles. We use brackets here, (...), to make 
explicit that this phcnomcnological expression deals with averages over macroscopically large numbers of particles, 
and to indicate that the particle concentration and flux can be meaningfully represented as smooth functions of space 
and time in macroscopic systems. Fick's First Law is the basis for Fick's Second Law, also known as the diffusion 
equation, 

m D m. (2) 

dt dx* ■ ( ' 

These equations have been extensively verified in bulk gases and solutions with macroscopically large number of 
particles Q. 

Our particular interest here is in the "small numbers" limit of Fick's Law, where there are only a few particles 
in the system, with special reference to the fluctuations that attend diffusive dynamics. Small particle numbers and 
their fluctuations are important in nanotechnology; inside biological cells, where the typical copy number of any 
given type of protein is often less than a few thousand 0|; and in single- molecule studies of ion channels, molecular 
motors, and in laser trap experiments 0,0, El- Fick's Law describes averages over a macroscopic number of particles; 
it does not describe small- number fiuctuational quantities, such as (J 2 ) — (J) 2 , or any other aspect of the flux 
distribution function. One of our motivations for undertaking this work is growing theoretical and experimental 
interest in nonequilibrium dynamics which centers on the distributions of microtrajectories available to such systems. 
We reasoned that a first step in examining the distribution of microtrajectories in nonequilibrium systems would be 
to revisit well-established problems such as diffusion, but with an eye to explicitly measuring (and calculating) the 
entire distribution of microscopic trajectories. 

Docs Fick's Law hold in the limit of small numbers of particles? And, are there violations? That is, if Fick's Law 
predicts flow to the right, due to a concentration gradient sloping downward towards the right, does it ever happen 
that particles flow instead to the left? Such situations have been called "Second-Law violations" 0,01; or, in classical 
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thermal problems, they are expressed in terms of "Maxwell's Demon" Such fluctuations are, of course, not real 
violations of the Second Law, because the Second Law is only a statement about averages, not fluctuations |lOj |. 
In this article we refer instead to such trajectories that go "against the grain" as bad-actors. Our interest here is 
not just in average fluxes, but in the full flux distribution function. Traditionally, diffusion has been studied in the 
bulk, where the number of particles is large. Only more recently has it become possible to perform experiments 
on small-numbers diffusion and on dynamical distribution functions, which explicitly emphasize the character of the 
microscopic trajectories that describe the dynamics of noncquilibrium systems, based on advances in nanotechnology, 
video microscopy and microfluidics. 

To predict dynamical distributions of diffusion rates, we can use either classical random -flig ht theory or a recent 
maximum-entropy-like approach ^lj, called maximum caliber, based on work of ET Jaynes |l2j . In short, if particles 
are independent, diffusing in one dimension, and if their jump rates are stationary in time, the distribution of particle 
fluxes, P(J), at time t along an x-axis from one bin at x having N\ particles, to an adjacent bin at x + Ax having 
N2 particles, should follow the binomial distribution, or approximately a Gaussian function 

p (j) _ ■ wnp /' 
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where A7V = Ni — N2, N = Ni + N2, (A J 2 ) is the variance in the flux J, and q = 1 — p, with p being the probability 
that a particle jumps in the time interval At. 

Various moments of the distribution function are readily obtained from this model. First, the model predicts that 
the average net number of particles, J, that jump per unit time at time t is [TH ] 

(J) = (:n-j2) = -j^AN, (4) 

where ji is the flux from the bin 1 at a; to bin 2 at x + Ax, and J2 is from bin 2 to 1. This proportionality of the 
average flux (J) to AN, simply predicts Fick's Law, where the diffusion coefficient D is related to p by D — pAx 2 /At, 
and where Ax is the bin size and At is the unit time step. 

For the flux fluctuations, i.e., the second moment, the model predicts 

(AJ 2 } = ((J-(J)f} = ^^-N, (5) 

where N = N% +N2 is the total number of particles associated with the two bins of interest. Hence, the key prediction 
here is that the flux fluctuations are proportional to the total particle number, N. 

We are also interested in the number of bad-actors, i.e., the number of trajectories that would lead to particle flows 
up a concentration gradient, rather than down it. This quantity can be derived from the flux distribution ^lj as 

^badactors = r I 1 ~ erf ' 
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(6) 



where the approximation holds for small values of (J) / yj ((AJ) 2 ). In the expression above, the next higher term 
(the cubic term), is an order of magnitude smaller than the linear term for the values of (J)/y ((A J) 2 ) used in our 
experiments (see Figure EJ. 



A. A "Flux" Fluctuation Theorem 



Recently, a useful description of noncquilibrium dynamics has involved fluctuation theorems. Fluctuation theorems 
characterize the extent to which the system deviates from its dominant flow behavior 0, 0, [TM . In the 
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FIG. 1: The microfluidics experiment. Colloids corralled on one side of a gate begin to diffuse at time t = by opening the 
gate, (a) Schematic of the microfluidic chip (see text for details), (b) The geometry of the microfluidic chamber (not drawn to 
scale) . 



diffusive dynamics case of interest here, if the number of particles Ni in bin 1 is greater than the number of particles 
N2 in bin 2, then particles, on average, will flow from 1 to 2. Fluctuation theorems describe the amount of reverse 
flow. Ours is a flux fluctuation theorem, i.e., it is expressed in terms of the quantity P(J)/P(—J), where J is the 
flux. This differs from fluctuation theorems @, El expressed in terms of entropies, P(AS)/P(—AS) and from the 
work theorems of Jarzynski and Crooks, which are expressed in terms of the work w as P[w/(kT)]/P[—w/(kT)] [l7j . 
In our approach, the ratio P(J)/P(—J), obtained from Eq. [31 gives the ratio of probabilities of fluxes in the forward 
and backward directions [Tl| . 

1 P ( J ) 2(J) 

Thus, the quantity ln[P(J)/P(— J)] is predicted to be proportional to the normalized flux (J)/((AJ) 2 ) x J. In situa- 
tions having large flux, the back-flow becomes exponentially negligible. We subjected these predictions to experimental 
tests. 



II. THE MICROFLUIDICS EXPERIMENTS 



To study the dynamical distributions in diffusion, we devised a microfluidics experiment. Using the techniques of 
soft lithography, chip fabrication [l8j and the Sylgard 184 Silicone elastomer kit (Dow Corning Corporation), we made 
a microfluidics chamber having approximate dimensions 400/im by 100/im, partitioned into two regions (see Fig. ^t). 
The cross-section of this chamber is a segment of a circular disc, with a maximum depth of 10/im (see Fig.^i). The 
chamber is filled on one side with a solution containing about 200 colloidal, green fluorescent polystyrene particles 
0.29/im in diameter (Duke Scientific. Cat. No. G300) (sec Fig.^). The beads are at an optimized concentration so 
that the interactions are negligible while at the same time permitting sufficient statistics over a wider range of 
AN and N. 

At time t = 0, we open a microfluidic gate, allowing particles to diffuse from one side to the other, and we begin 
taking periodic snapshots under an Olympus 1X71 inverted microscope. (We performed the same experiment under 
equilibrium conditions where the initial concentration was uniform across the whole chamber (results not shown, 
see H3))- We take 3 snapshots of the beads in the chamber every time interval of At = 10 seconds for 6 hours. 
Since there is a possibility that within one snapshot some particles are temporarily overlapping and/or are out of 
focus, the 3 snapshots arc used to minimize that error. The snapshots are taken using fluorescence microscopy with 
a SONY DFW-V500 camera. (During the time when no snapshots are taken, a shutter prevents the experimental 
chamber from being exposed to the incident light, to prevent photobleaching and heating the chamber.) We then 
determine the particle positions at each snapshot using a computerized centroid tracking algorithm |25j . The time- 
dependent particle density is determined by dividing the chamber into a number of equal-sized bins of value Ax each 
along the longest dimension of 400/im and by computing the number of particles in each bin as a function of time. 
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Although the microfhiidic chamber is three-dimensional, it can be shown that, in the case of weak particle-particle 
and particle-wall interactions, the problem can be collapsed to a one-dimension diffusion problem. Therefore, we bin 
only along the x-axis, the direction of the concentration gradient. The choice of the bin size will affect the statistics 
for each combination of N± and N2 as well as the range of N and AJV themselves. If the bin size is too big, then 
there will not be enough statistics and, in addition, the range of values for ./V will not include small numbers (since 
it will be rare to have 1 or 2 particles in a single bin). On the other hand, if the bin size is too small, we may 
not have a sufficient range of values for N and AiV (since it will be rare to have more than a couple of particles 
in the small bin). Also, in this case, there will be an increased probability for a particle to have multiple jumps 
across bins within the time interval At. Therefore, the optimal choice of the bin size was made based on the bead's 
expected mean excursion within the time-interval At, which is \J1DAt. This is the only relevant microscopic length 
scale. Here, D is the diffusion coefficient for an individual bead given by the Stokes's formula [2^. For a bead of 
0.29^m in diameter suspended in water at room temperature, the Stokes's formula gives a diffusion coefficient D of 
approximately 1.5/im 2 /second. This value, within experimental error, is equal to the one we obtain by fitting our 
data of the concentration profile at different times to the one dimensional diffusion equation using D as our fitting 
parameter (i.e. D = 1.3 ± 0.27^m 2 /second). This gives a bin size of Ax r* By observing all the consecutive bin 

pairs for all the frames taken we were able to obtain, on average, about 5000 points for each combination of N% and 
N2- Given the bead concentration in the microfluidics channel, the N\ and N2 ranged from through 6. The choice 
of bin size determines the value of the jump probability p, as discussed in |2lj . 

We can find the flux at a plane i at a specific time interval from the computed particle distribution statistics as 
a function of position x and time t mentioned above. Since the microfluidic chamber is isolated, the total number 
of particles stays the same from one frame to the next. As a result of this conservation in particle number, the flux 
at plane i + 1 (Jj+i), i.e., the plane that separates bins i and i + 1, can be easily evaluated by using the continuity 
equation, 

Ni(t + At) = Ni(t) + (Ji(t)-J i+1 (t))At (8) 

where iVj is the number of particles in bin i. Since the microfluidic chamber is isolated, from our boundary conditions 
the flux Jo (flux at x = 0) is zero at all times. Combined with Eq.El we obtain Ji(i). Thus, from the analysis of these 
images, we obtain complete sets of the values of {iVj(i)} and {Ji(t)} in all the bins and at all times of observation. 
Then, for each pair of consecutive bins with specific values of N\ and N%, we construct the histogram of J values. 
Upon normalization the histogram becomes the flux probability distribution, P(J). 

III. RESULTS 

A. The Flux Distribution Function is Gaussian. 

Figure [3 shows our observed particle flux distribution function at the optimized concentration. All the data falls on 
a single master curve where (J) and ((AJ) 2 ) have been calculated separately from each combination of N\ and N2. 
The quadratic form observed on this log plot shows that the distribution function is Gaussian. The theory predicts 
that: (i) the coefficient of the square term should be — 1, (ii) the coefficient of the linear term should be zero, and 
(iii) the constant term should be ln(A Jbi n /\/27r) ~ —0.9, where AJbm is the bin-size used to obtain the histogram 
and is equal to 0.1 second -1 , i.e., 1 particle per unit time. Consistent with these predictions, the coefficient observed 
for the square term is —0.98, for the linear term is —0.0018, and for the constant term is —0.94. The coefficient of 
determination for the quadratic fit is R 2 = 0.98. Next, we analyze the bad actors - the backward flows - in two 
different ways. 

B. The Bad-Actor Trajectory Counts are Well Predicted by the Model. 



Equation[f)]predicts that, for small values of (J) / y/ ((A J) 2 ), the fraction of bad-actors should be linearly proportional 
to (J)/y/ ((AJ) 2 ). In good agreement, Figure ^confirms this linearity and gives the predicted intercept of 0.5. This 
means that as the system approaches equilibrium (i.e. (J) w 0), about half the trajectories involve flow down the 
vanishingly small gradient and half the trajectories involve flow up that small gradient. In the linear regime, the 
best fit line shows the slope to be 0.37, which agrees well with the expected value of w 0.4 from Eq. |SJ The 

coefficient of determination for the linear fit is R 2 = 0.99. Another key feature of this graph is that when the system 



5 




-1 o 1 

(J - <J>)/(2<AJ 2 >) a5 



FIG. 2: The flux distribution function. \ ln(((AJ) 2 )) + ln(P(J)) is plotted against (J - (J))/^/2{(AJ) 2 ), based on the form 
indicated by Equation [3] Circles indicate experimental points, the line shows a quadratic fit to the data. The coefficient of 
determination (R 2 ) for the fit is also reported. This demonstrates that the distribution function is Gaussian, and we find that 
the coefficients are well predicted by Equation 
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FIG. 3: The fraction of trajectories that are bad-actors vs. the deviation from equilibrium as characterized by the normalized 
mean flux, (J) / \J ((A J) 2 }. Experimental data is shown in squares, while the solid line represents the fit to the data. The 
coefficient of determination (R 2 ) for the fit is also reported. The slope and intercept agree well with the model. 



is farther away from equilibrium (as implied by a larger mean flux), the bad-actor fraction reduces. What this means 
is that more of the microtrajectories available to the system in this case are potent to change the current state of the 
system. Indeed, in our earlier paper we characterized this idea quantitatively in the form of the potency \H\ . 



C. The Experiments Confirm the Flux Fluctuation Theorem. 



Figure 0] shows hx(P(J)/P(—J)) vs. the flux, normalized by (J) J/((AJ) 2 ), to account for different averages and 
variances of the flux distribution. This rescaling leads to a linear master curve as predicted by Eq. Experiments 
show the slope to be 2.0, in perfect agreement with the predicted slope of 2 from Eq.0 In this figure, there are a few 
aberrant data points that we are unable to explain. However, the coefficient of determination is still close to unity, 
R 2 = 0.77. Though ultimately the fluctuation theorem is a reflection of the nature of the flux distribution function, 
such theorems are a compact way to quantitatively illustrate the significance of bad-actor microtrajectories. 
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FIG. 4: The flux fluctuation theorem. The plot shows ln[P(J)/P(—J)] vs. (J) J/{(AJ) 2 ) for different values of (J) and <(AJ) 2 } 
arising due to different combinations of Ni and N2- Experimental data is shown in circles, while the solid line represents the fit 
to the data. The coefficient of determination (R 2 ) for the fit is also reported. The slope and intercept agree with the prediction 
of Equation [7| 




AN 

FIG. 5: Experimental support for Fick's Law, even down to few-particle gradients. The average flux, (J), is shown as a function 
of AN, the gradient in the particle number between two neighboring bins. Experimental data is shown in circles, while the 
solid line represents the fit to the data. The coefficient of determination (R 2 ) for the fit is also reported. The error bars shown 
are the variances due to the different combinations of Ni and N2 resulting in the same AN. The slope and intercept are in 
agreement with the expected theoretical values, based on Equation 0] 

D. Fick's Law Holds Even in the Small-Numbers Limit. 

We compare the average flux between two neighboring bins, (J), with the difference in particle numbers, AN = 
Ni — N2. This data is compiled from all the values of N\ and A^ that provide a given AN. Fig. [3] shows that (J) 
depends linearly on the particle number gradient AN, even down to "gradients" of zero or one particle, indicating 
that Fick's Law holds in the small-numbers limit. The slope of the graph (0.03/second) also gives us a value of the 
jump rate, p = 0.3, which is in good agreement with the theoretical estimate of 0.33 made in terms of the bin-size and 
the diffusion coefficient |2l|. As expected, the intercept is close to (see Equation QJ. The linear fit has a coefficient 
of determination R 2 = 0.99. 
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FIG. 6: The second cumulant, (AJ 2 ) = (J 2 ) — (J) 2 , vs. the total number of particles, N. The second moment of particle 
flux is proportional to the sum of particle numbers in the two bins. Experimental data is shown in circles, while the solid 
line represents the fit to the data. The coefficient of determination (R 2 ) for the fit is also reported. The error bars show the 
variances due to the different combinations of Ni and N2 that result in the same N. The slope and intercept are well predicted 
by Equation 



E. The Second Moment of Particle Flux is Proportional to the Sum of Particle Numbers. 

Equation predicts that the second moment of the flux should be proportional to the sum of particle numbers in 
the two bins, N = Ni + iVjj. Figure O confirms this dependence of (AJ 2 ) = ((J — (J)) 2 ) on N. For the optimized 
particle concentrations, the slope (0.0022/second 2 ) is equal to the expected slope of 0.0022/second 2 for the value of 
p = 0.33. The coefficient of determination is R 2 = 0.96. At higher particle concentrations (data not shown), however, 
not surprisingly, systematic errors begin to appear and the slope deviation is quite high compared to the expected 
value. We performed Brownian dynamics simulations that show the likely cause of these concentration-dependent 
errors is non-conservation of bin counts, from particles that either overlap or go out of focus in one snapshot and into 
focus in the next (see previous section) . 

IV. CONCLUSIONS 

In a microfiuidics experiment, we have determined the distribution of particle fluxes in few-particle diffusion. 
Random-flight and maximum-caliber models predict a gaussian distribution of fluxes. Moreover, with a single pa- 
rameter p, which is essentially the diffusion constant, we find agreement of the theory with several experimental 
properties which are usually not examined in diffusion, and are in direct analogy with quantities of recent interest in 
other nonequilibrium experiments [1,0, 0,0]- First, we find that Fick's law - the proportionality of average flux to 
the gradient of average concentration - holds even down to concentration gradients as small as a single particle. Ex- 
periments also confirm that the variance in the flux is proportional to the total number of particles, (J 2 ) oc N\ + N%, 
with correct slopes within experimental errors. In addition, we describe a new"flux fluctuation theorem", that is 
found to be consistent with the data in predicting an exponentially diminishing number of variant trajectories, as 
a function of the deviation from equilibrium. The model predicts the backwards flows, the bad actors, which are 
relatively infrequent situations in which particles flow up — rather than down — their concentration gradients and 
shows that this subset of the overall repertoire of microscopic trajectories can be characterized quantitatively. 
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